#!/usr/bin/env python
"""
Vizualize mesh partitionsi using gnuplot or vtk

GNUPLOT:
Build a gnuplot data file and the corresponding gnuplot script (need gnuplot installed)

VTK:
Render the scene on screen (need vtk installed)

INPUT
-n triangle
Two input files describing the mesh are needed:
.node contains the nodes
.edge contains the edges
These two files are required

-n gts
One input file only is required to describe the mesh
.gts contains the nodes coordinate, edge and faces (useless here) description
This file is required

-n gmsh
One input file is required to describe the mesh
.msh contains the nodes coordinates and elements description
This file is required

-n graph
2 input files
.node with the node coordinates (same as metis/triangle)
.graph: connectivity graph description

-i triangle
.part in the triangle format (see below \"FORMAT\"). Contains the domain assignment for each node

-i metis/pagrid
.part in the metis/pagrid format (see below \"FORMAT\")
The part files are not required (if present, the partition will be drawn)


OUTPUT
Render: -r
-r gnuplot
Write two different files:
- .gpd: contains the data to draw
- .gps: contains the gnuplot script that read the gnuplot data text
By default, these files are deleted (option -t), kept with -p
2 (3) output files:
--pdf: pdf output (uses ps2pdf, not direct pdf)
--png: png output (for internet display, uses imageMagick's convert option of .ps, no direct output (gnuplot'gd is ugly))
--svg: W3C svg output. For internet display. Currently only working good with Opera, slight differences with Omniweb, not good at all for Firefox and Safari


-r vtk. Uses Visual ToolKit
Warning: Long time to load vtk python library. Long time to render if number of vertices >1000. Created files can be very big (>50 megas for 1740 vertices)
--interactive: Interactive use. Possibility of saving with a simple printScreen (DEFAULT)
--pdf: 2D pdf output. Try hard the pdf 3D (using meshlab, U3D, movie15's pdflatex) but not working for now. Seems to be viewable only with Adobe reader for now (no preview.app, or xpdf, ...)
--rasterpdf: raster pdf. Smaller than normal pdf, but loose quality. Means to be faster rendering.
--png: png image file format. Internet use.
--vrml: W3C virtual world standard. Use with internet and correct plug-in (freewrl, octoga, cosmo player, on windows flux, ...). X3D not possible for now.
--obj: OBJ Wavefront file format (2 files: .obj referencing a .mtl). Means to be an exchange format between renderers. 


USAGE:

showmesh2.py [options] mesh file(s) [partition file]
All the options have to be BEFORE the 3 input files

OPTIONS:

gnuplot options start with --gpl
vtk options start with --vtk

-h, --help Print this help
-t Delete temp files (.gpd and .gps). Dont use this option for now, since the automatic plotting does NOT work
-p Do NOT delete temp files
--nopart Dont draw partitions, even if a .part file is present (DEFAULT: partitions drawn if .part present, not drawn if not)
--paper Use templates more adapted for a paper thesis
--gpl_out_data arg .gpd file name
--gpl_out_script arg .gps file name
-o, --output arg prefix name of the output file (suffix can be .eps if gnuplot renderer, or .pdf/.png if vtk renderer)
-d, --debug Print some debug infos
-i, --in_format arg Input format (for part file in particulat): arg can be triangle or metis (DEFAULT=Metis)
-r, --render arg Output renderer: arg can be gnuplot or vtk (DEFAULT=gnuplot)
--rsphere arg Specifies radius of the nodes (for gnuplot and vtk renderers) (DEFAULT: 0.6 for gnuplot AND VTK)
--rtube arg Specifies radius of the edges (DEFAUTL: vtk_rsphere/10 for VTK, 0.6 for gnuplot)
--pdf, --png, --svg, --interactive (VTK), --rasterpdf (VTK), --vrml(VTK), --obj(VTK)  Format of the output file (DEFAULT: vtk_nooutput). vtk_inter means interactive vizualization
-l, --label: print node number. Start from 1 (metis standard). Size of the text corresponds to the size of the sphere, so do not use for big graphs, but rather for debugging purposes (implied in -d (debug) option)


FILES FORMAT

-i metis/triangle
.node format
# #_of_nodes dimension (2 or 3)
# node_id x y (z)
9 2 #number of nodes, dimension
1 0 0 0 # node_number, x, y, z (if dim=3)
2 1 0 0
3 2 0 0
4 0 1 0
5 1 1 0
6 2 1 0
7 0 2 0
8 1 2 0
9 2 2 0

.edge format
#Edges
5 0 #Number of edges, 0
1 1 5 #A number, node1, node2
2 8 9
3 5 8
4 9 6
5 2 3


.part format (Triangle)
# hophophop
9 4 #Number of nodes, number of domains
1 1 # Node, domain
2 1
3 4
4 1
5 1
6 3
7 2
8 2
9 2

.part format (Metis)
#hiphip
1 # domain #. ith line represents node i
2
2
1
4
3
4
4
2

EXAMPLE:

vtk:
  ./showmesh.py -n triangle -i metis -r vtk --rsphere 0.6 -o test1 --vrml ../data/testix.1.node ../data/testix.1.edge ../data/parts_metis/testix.1.mesh.ngraph.part.7

Basile Clout, September 2007-February 2008
"""


import os
import sys
import string
import getopt
import pprint
from subprocess import *
from datetime import *
from string import Template

def pnode(fnode):
    """ Parse the .node file
    and returns the correspondind
    list of the nodes
    """

    nodes = []
    count = 0
    dim = 0
    nbnodes = 0

    for line in fnode:

        if line[0]=='#':
            continue
        else:
            parts = line.split()

            if count == 0:
                nbnodes = int(parts[0])
                dim = int(parts[1])
                count += 1
                if dim not in [2,3]:
                    print "Error dim in .node file!"

            else:
                if dim == 3:
                    nodes.append(tuple([float(x) for x in parts[1:]]))
                else:
                    nodes.append(tuple([float(x) for x in parts[1:3]]))

    if (len(nodes) != nbnodes):
        print "Error: actual number of nodes does not corresponds " + str(len(nodes)), nbnodes
        sys.exit(2)
                                 
    return nodes, dim


def gts_node(fnode):
    """ Parse the .node file in a gts way
    Format: First line: #node(nv) #edges(ne) #elements(nf)
    nv next lines: 3D coordinates of the node
    ne followng lines: 2 nodes indices cooresponding to the vertices of the edge
    nf following lines: 3 nodes indices corresponding to the vertices of the triangular face (dont care here)
    """
    
    nodes=[]
    edges=[]
    cnt_line = 0
    nv = 0
    ne = 0
    nf = 0

    
    for line in fnode:

        if line[0]=='#':
            continue

        parts = line.split()

        if cnt_line == 0:
            (nv, ne, nf) = (int(x) for x in parts)

        else:
            if (cnt_line < nv+1):
                nodes.append(tuple([float(x) for x in parts]))
            elif (cnt_line < ne+nv+1):
                edges.append(tuple([int(x) for x in parts]))
            else:
                break

        cnt_line+=1
        
    if len(edges) != ne:
        print "ERROR: Actual number of edges does not correspond "+str(len(edges)), str(ne)
        sys.exit(2)
    if len(nodes) != nv:
        print "ERROR: Actual number of nodes does not correspond "+str(len(nodes)), str(nv)
        sys.exit(2)

    return nodes, edges


def gmsh_node(fnode):
    """ Parse the .msh file in a gmsh way
    Format: Skip everything (header) until reaching the '$Nodes' keyword
    Next line = number of nodes
    Next lines = node# (starting from 1) x y z
    Reach $EndNodes
    Reach $Elements
    Next line = number of elements
    Net lines = element# geometrical_type (15=1-node point, 1=2-nodes line, 2=3-nodes triangle in the simple first order case) tags node1 node2 node3 ...#
    Reach $EndElements
    """

    nodes=[]
    edges=[]
    nbnodes=0
    nbele=0
    mod = 0                            # 0=nothing, 1=nodes, 2=elements
    cnt = 0

    for line in fnode:

        parts = line.split()
        keyword = parts[0]

        if keyword == "$Nodes":
            mod = 1
            cnt = 1 
            continue
        elif keyword == "$Elements":
            mod = 2
            cnt = 1
            continue
        elif keyword in ("$EndNodes", "$EndElements"):
            mod=0
            continue
        else:
            if mod == 0:
                continue
            else:
                if cnt == 1:
                    if mod == 1:
                        nbnodes = int(keyword)
                    elif mod == 2:
                        nbele = int(keyword)
                    else:
                        print "Error with mod"
                        sys.exit(2)
                    cnt=0
                else:
                    if mod == 1:
                        nodes.append(tuple([float(x) for x in parts[1:]]))
                    elif mod == 2:
                        ele_type = int(parts[1])
#                        if ele_type == 1:
#                            edges.append(tuple([int(x) for x in parts[6:]]))
                        if ele_type == 2:
                            (n1, n2, n3) = (int(parts[6]), int(parts[7]), int(parts[8]))
                            if (n1, n2) not in edges: 
                                edges.append((n1, n2))
                            if (n2, n3) not in edges:
                                edges.append((n2, n3))
                            if (n1, n3) not in edges:
                                edges.append((n1,n3))
        
    if len(nodes) != nbnodes:
        print "ERROR: Actual number of nodes does not correspond "+str(len(nodes)), str(nbnodes)
        sys.exit(2)

    return nodes, edges
        

        
def pedges(fedge):
    """ Parse the .edge file and returns the corresponding
    list of edges
    """

    edges = []
    count = 0
    nbedges = 0
    boundary = 0
    for line in fedge:

        if line[0]=='#':
            continue
        else:
            parts = line.split()

            if count == 0:
                nbedges = int(parts[0])
                boundary = int(parts[1])
                count += 1
                if boundary not in (0,1):
                    print "Error boundary in .edge file!"

            else:
                if boundary == 1:
                    edges.append(tuple([int(x) for x in parts[1:]]))
                else:
                    edges.append(tuple([int(x) for x in parts[1:3]]))

    if (len(edges) != nbedges):
        print "Error: actual number of edges does not corresponds " + str(len(edges)), str(nbedges)
        sys.exit(2)

    return edges


def gedges(fedge):
    """ Parse a graph description of a mesh (input for partitioners
    for example) and return the corresponding list of edges
    """

    edges=[]
    cnt=0
    nbedges=0
    nbelements=0

    for line in fedge:

        parts=line.split()
        if cnt==0:
            (nbelements,nbedges)=(int(parts[0]), int(parts[1]))
        else:
            for i in range(0,len(parts)):
                if ((cnt, int(parts[i])) not in edges) and ((int(parts[i]), cnt) not in edges):
                    edges.append((cnt, int(parts[i])))
        cnt+=1

#    if len(edges)!=nbedges:
#        print "ERROR: Actual number of links does not correspond "+str(len(edges)), str(nbedges)
#        sys.exit(2)
    return edges

    

def pparts(fparts, fl_format):
    """ Build a list of the partitions:
    [(node, part), (node, part), ...]
    """

    nbnodes = 0
    nbparts = 0
    count = 0
    mapping = []
    domains = []
    i = 0
    p = 0
    dom = 0

    for line in fparts:

        if line[0]=='#':
            continue
        
        else:
            parts = line.split()

            if fl_format == 0: # Check the first line if in triangle format
            
                if count == 0 :
                    nbnodes = int(parts[0])
                    nbparts = int(parts[1])
                    count += 1
                    
                else:
                    if int(parts[1]) not in domains:
                        domains.append(int(parts[1]))
                        
                    mapping.append(tuple([int(x) for x in parts]))
                        
            elif fl_format == 1:

                dom = int(parts[0]) + 1
                if dom not in domains:
                    domains.append(dom)
                    p = p+1

                i = i+1
                mapping.append((i, dom))

    if(fl_format==1):
        nbnodes = i;
        nbparts = p;
                    
    if (len(mapping) != nbnodes ):
        print "Error: actual number of nodes does not corresponds " + str(len(mapping)), nbnodes
        sys.exit(2)
    if (len(domains) != nbparts):
        print "Error: actual number of domains does not corresponds " + str(len(domains)), nbparts
        sys.exit(2)

    domains.sort()
    return mapping, domains


def edgeType(edges, mapping):
    """ Return true if the edge is a cutting edge
    """
    cutedges = []
    oedges = []
    for edge in edges:
        if mapping[edge[0]-1][1] != mapping[edge[1]-1][1]:
            cutedges.append(edge)
        else:
            oedges.append(edge)

    return cutedges, oedges
    

def pcoord(coord, dim):
    """ Returns a string containing the formatting of
    the coordinates from coord, given the dimension (2 or 3)
    """

    if dim == 3:
        return str(coord[0]) + " " + str(coord[1]) + " " + str(coord[2]) +"\n"
    else:
        return str(coord[0]) + " " + str(coord[1]) + "\n"



def build_gpd(nodes, dim, edges, domains, mappings, out):
    """ Build the gnuplot data file
    Returns the string to print in this file
    """

    intro = "# Gnuplot data for mesh vizualization\n"
    intro += "# Generated by " + sys.argv[0] + " on " + str(datetime.now()) + " from the files " + str(sys.argv[-3:]) + "\n"
    intro += "# The Nodes paragraph contains the nodes, the Edges one the ... edges description\n\n"
    
    intro += "# Nodes: id x y z\n"
    out.write(intro)

    for i in domains:
        out.write("# Domain " + str(i) + "\n")
        for (node, domain) in mappings:
            if i==domain:
                out.write(pcoord(nodes[node-1], dim))

        out.write("\n\n")
        
    
    out.write("# Edges: node1 node2\n")

    # Cutedges list contains edges between two nodes within two different domains
    cutedges, oedges = edgeType(edges, mappings) 
    #print "cutedges:"
    #pprint.pprint(cutedges)
    #print "other edges:"
    #pprint.pprint(oedges)
    
    intro += "# No cutting-edges\n"
    for edge in oedges:
        idnode1 = nodes[edge[0]-1]
        idnode2 = nodes[edge[1]-1]
        out.write(pcoord(idnode1, dim))
        out.write(pcoord(idnode2, dim))
        out.write("\n")

    out.write("\n")
    out.write("# Cutting edges\n")

    if cutedges != []:
        for edge in cutedges:
            idnode1 = nodes[edge[0]-1]
            idnode2 = nodes[edge[1]-1]
            out.write(pcoord(idnode1, dim))
            out.write(pcoord(idnode2, dim))
            out.write("\n" )
            

    return cutedges

def build_gps(nodes, dim, edges, domains, mappings, filename, out_gpd, out, fl_paper, ptsize, cutedges, rtube, fl_type, fl_legend, names):
    """ Build the gnuplot script file
    Returns the string to print in this file
    """

    intro = "# Gnuplot script for mesh vizualization\n"
    intro += "# Generated by " + sys.argv[0] + " on " + str(datetime.now()) + " from the files " + str(sys.argv[-3:]) + "\n"
    if fl_type in (1,3):
        intro += "set term postscript color\n"
        intro += "set output \"%s.ps\"\n" % filename
    elif fl_type == 6:
        intro += "set term svg\n"
        intro += "set output \"%s.svg\"\n" % filename
    else:
        print "Bad fl_type"
        sys.exit(1)
    intro += "set size ratio -1\n"      # Orthogonal plot
    if fl_legend == 1:
        intro += "set key right bottom\n"
    else:
        intro += "set key off\n"
        
    if len(domains)==1:
        title = "Graph %s" % names[0]
    else:
        title = "Partition %s (graph %s)" % (names[1], names[0])
        
    intro += "set title \'"+title+"\'\n"
    intro += "set pointsize %d\n" % ptsize
    intro += "unset border\n"
    intro += "unset xtics\n"
    intro += "unset ytics\n"
    xmin, xmax, ymin, ymax = window(nodes)
    intro += "set xrange [%d:%d]\n" % (xmin, xmax)
    intro += "set yrange [%d:%d]\n" % (ymin, ymax)
    if (fl_paper==1):
        intro+="set palette rgbformulae 0,0,0\n"
    
    out.write(intro)

    out.write("\n# Plot\n")
    if dim==3:
        plot = "splot "
        myrange = "1:2:3"
    else:
        plot = "plot "
        myrange = "1:2"
        
    out.write(plot)

    tnodes = Template("\'" + out_gpd +"\' index $index using $urange title '$title' with points pt $pt ps $ray lc $color$end")
    tedges= Template("\'" + out_gpd +"\' index $index using $urange title '$title' with l lw $lwsize lt $lt lc $color$end")

        
    colors_ps=[0, 1, 8, 6, 4, 3, 5, 2]
    cc_ps = 9
    co_ps = 1
    colors_svg=[-1, 1, 7, 9, 12, 3, 4, 2]
    cc_svg = 0
    co_svg = 1

    if fl_type in (1,3):
        mycolors = colors_ps
        cc = cc_ps                      # Grey color for postscript
        co = co_ps                      # Red color for postscript
    elif fl_type == 6:
        mycolors = colors_svg
        cc = cc_svg
        co = co_svg
        
    e = len(domains)                    # First draw domains.
    coef_cc = 5                         # ratio thickness cutedge/otheredge

    # Print edges
    if fl_paper==0:                     # Normal output
        out.write(tedges.substitute(index = e, color = cc, lt = 1, urange = myrange, lwsize = rtube, title = "domain edges", end = ",\\\n"))
        if cutedges != 0:
            out.write(tedges.substitute(index = e+1, color = co, lt = 1, urange = myrange, lwsize = coef_cc*rtube, title = "cutting edges", end = ",\\\n"))
    else:
        out.write(tedges.substitute(index = e, color = cc, lt = 2, urange = myrange, lwsize = rtube, title = "domain edges", end = ",\\\n"))
        if cutedges != 0:
            out.write(tedges.substitute(index = e+1, color = co, lt = 3, urange = myrange, lwsize = coef_cc*rtube, title = "cut edges", end = ",\\\n"))
    
    # Print nodes
    if len(domains) == 1:
        out.write(tnodes.substitute(index = 0, color = 0, pt = 7, urange = myrange, title = "nodes", ray=ptsize, end = ""))
    else:
        for i in domains:
            if fl_paper == 0:
                if i==len(domains):
                    out.write(tnodes.substitute(index = i-1, color = mycolors[i-1], pt = 7, urange = myrange, ray = ptsize, title = "processor " + str(i-1), end = ""))
                else:
                    out.write(tnodes.substitute(index = i-1, color = mycolors[i-1], pt = 7, urange = myrange, ray = ptsize, title = "processor " + str(i-1), end = ",\\\n"))
            else:
                if i==len(domains):
                    out.write(tnodes.substitute(index = i-1, color = 1, pt = 7, urange = myrange, ray = ptsize, title = "processor " + str(i), end = ""))
                else:
                    out.write(tnodes.substitute(index = i-1, color = 1, pt = 7, urange = myrange, ray = ptsize, title = "processor " + str(i), end = ",\\\n"))
            
#    out.write("\n\n print \"" + filename + " created!\"\n");


def vtk_label_node(i, node, rsphere):
    """ Insert the 3D label corresponding to node
    """

    import vtk
    
    if len(node)==2:
        attachment = [node[0], node[1], 0]
    else:
        attachment = node
        
    caption1 = vtk.vtkCaptionActor2D()
    caption1.GetPositionCoordinate().SetCoordinateSystemToWorld()
    caption1.GetPositionCoordinate().SetValue(-rsphere/2, -rsphere/2, -rsphere/2)
    caption1.SetWidth(0.1)
    caption1.SetHeight(rsphere/100)
    caption1.BorderOff()
    caption1.ThreeDimensionalLeaderOn()
    caption1.SetAttachmentPoint(attachment)
    caption1.SetCaption(str(i+1))
    caption1.GetCaptionTextProperty().SetFontSize(10)

    return caption1


def build_vtk(_nodes, _dim, _edges, _domains, _mappings, _rsphere, _rtube, filename, fl_file, fl_legend, _label):
    """ Build vtk rendering from standard objects
    """

    import vtk

    color = "range"                     # set to "fix" if use the following color's table
    
    xmin, xmax, ymin, ymax = window(_nodes)
    wx, wy = 1000, 700                   # Rendering window size (and size of the printed graphics)
    labels = []
    scale=[_rsphere*2, _rsphere*2, _rsphere*2]

    colors = [   
               [0.5, 0.5, 0.5, 1.0],    # Grey
               [0.2, 0.2, 1.0, 1.0],    # Deep Blue
               [0.5, 0.0, 0.3, 1.0],    # Deep red
               [0.0, 1.0, 0.0, 1.0],    # Green
               [1.0, 1.0, 0.0, 1.0],    # Yellow
               [1.0, 0.5, 0.0, 1.0],    # Orange
               [1.0, 0.0, 1.0, 1.0],    # Magenta
               [0.0, 0.6, 0.2, 1.0],    # Dark green
               [0.0, 0.0, 0.0, 1.0],    # Black
               [1.0, 1.0, 1.0, 1.0],    # White
               [1.0, 0.0, 0.0, 1.0],    # Red
             ]


  
    lut = vtk.vtkLookupTable()
    if color == "range":
        lut.SetNumberOfColors(512)
        lut.SetHueRange(0.08,0.92)
#        lut.SetSaturationRange(0.5, 1)
#        lut.SetValueRange(0.6, 1)
        lut.Build()

    elif color == "fix":
        lut.SetNumberOfColors(len(_domains))
        lut.Build()
        for i in range(0,len(_domains)):
            lut.SetTableValue(i, colors[i])
            
    lutEdges = vtk.vtkLookupTable()
    lutEdges.SetNumberOfColors(2)
    lutEdges.Build()
    lutEdges.SetTableValue(0, colors[0])  # Grey for domain edges
    lutEdges.SetTableValue(1, colors[10]) # Red for cut edges
    
    # Nodes (Points)
    nodes = vtk.vtkPoints()
    nodes.SetDataTypeToFloat()
    nodesColor = vtk.vtkIdTypeArray()
    nodesCell = vtk.vtkCellArray()
    i=0
    
    for node in _nodes:

        if(_dim==2):
            id = nodes.InsertPoint(i, node[0], node[1], 0.0)
        else:
            id = nodes.InsertPoint(i, node[0], node[1], node[2])
        
        nodesColor.InsertNextValue(_mappings[i][1])    # Color (domain) corresponding to this node

        nodesCell.InsertNextCell(1)
        nodesCell.InsertCellPoint(i)

        if _label==1:
            mylabel = vtk_label_node(i, node, _rsphere)
            labels.append(mylabel)

        i = i+1
        
    # Add nodes corresponding to the legend spheres
    if(fl_legend!=0):
        
        for part in _domains:
            id = nodes.InsertPoint(i, xmax+(xmax-xmin)/10, ymin+(part-1+2)*(ymax-ymin)/10, 0.0)
            nodesColor.InsertNextValue(part)
            nodesCell.InsertNextCell(1)
            nodesCell.InsertCellPoint(i)
            i = i+1
        
    nodesPoly = vtk.vtkPolyData()
    nodesPoly.SetPoints(nodes)
    nodesPoly.GetPointData().SetScalars(nodesColor)
   # nodesPoly.SetVerts(nodesCell)

    sphere = vtk.vtkSphereSource()
    sphere.SetThetaResolution(12)        # Depends on the power of the computer...
    sphere.SetPhiResolution(12)
    sphere.SetRadius(_rsphere)

    spheres = vtk.vtkGlyph3D()
    spheres.SetInput(nodesPoly)
    for i in _domains:
        spheres.SetSource(i-1, sphere.GetOutput())
    spheres.SetColorModeToColorByScalar()
    spheres.SetScaleModeToDataScalingOff()
#    spheres.SetIndexModeToScalar()
#    spheres.SetRange(1, len(_domains)*2)
#    spheres.SetRange(0, len(_colors))

    #Edges
    cutedgeColors = vtk.vtkIdTypeArray()
    cutedgeArray = vtk.vtkCellArray()
    cutedgePoints = vtk.vtkPoints()
    cutedgePoints.SetDataTypeToFloat()

    cutedges, oedges = edgeType(_edges, _mappings)
    
    for edge in cutedges:
        if _dim == 2 :
            pt1 = cutedgePoints.InsertNextPoint(_nodes[edge[0]-1][0], _nodes[edge[0]-1][1], 0.0)
            pt2 = cutedgePoints.InsertNextPoint(_nodes[edge[1]-1][0], _nodes[edge[1]-1][1], 0.0)
        else:
            pt1 = cutedgePoints.InsertNextPoint(_nodes[edge[0]-1])
            pt2 = cutedgePoints.InsertNextPoint(_nodes[edge[1]-1])

        cutedgeColors.InsertNextValue(12)

        cutedgeArray.InsertNextCell(2)
        cutedgeArray.InsertCellPoint(pt1)
        cutedgeArray.InsertCellPoint(pt2)


    if fl_legend != 0:
        pt1 = cutedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10, ymin, 0.0)
        pt2 = cutedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10+(xmax-xmin)/22, ymin, 0.0) 

        cutedgeColors.InsertNextValue(12)
        cutedgeArray.InsertNextCell(2)
        cutedgeArray.InsertCellPoint(pt1)
        cutedgeArray.InsertCellPoint(pt2)
        
    cutedgesPolyData = vtk.vtkPolyData()
    cutedgesPolyData.SetPoints(cutedgePoints)
    cutedgesPolyData.SetLines(cutedgeArray)
    cutedgesPolyData.GetCellData().SetScalars(cutedgeColors)


    # Other edges
    oedgeColors = vtk.vtkIdTypeArray()
    oedgeArray = vtk.vtkCellArray()
    oedgePoints = vtk.vtkPoints()
    oedgePoints.SetDataTypeToFloat()

    if fl_file != 3:
        oecolor = 0
    else:
        oecolor = 0
    
    for edge in oedges:
        if _dim == 2 :
            pt1 = oedgePoints.InsertNextPoint(_nodes[edge[0]-1][0], _nodes[edge[0]-1][1], 0.0)
            pt2 = oedgePoints.InsertNextPoint(_nodes[edge[1]-1][0], _nodes[edge[1]-1][1], 0.0)
        else:
            pt1 = oedgePoints.InsertNextPoint(_nodes[edge[0]-1])
            pt2 = oedgePoints.InsertNextPoint(_nodes[edge[1]-1])

        oedgeColors.InsertNextValue(oecolor)
        
        oedgeArray.InsertNextCell(2)
        oedgeArray.InsertCellPoint(pt1)
        oedgeArray.InsertCellPoint(pt2)

    if fl_legend != 0:
        pt1 = oedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10, ymin+(ymax-ymin)/10, 0.0)
        pt2 = oedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10+(xmax-xmin)/22, ymin+(ymax-ymin)/10, 0.0) 
        oedgeColors.InsertNextValue(oecolor)
        oedgeArray.InsertNextCell(2)
        oedgeArray.InsertCellPoint(pt1)
        oedgeArray.InsertCellPoint(pt2)
        
    oedgesPolyData = vtk.vtkPolyData()
    oedgesPolyData.SetPoints(oedgePoints)
    oedgesPolyData.SetLines(oedgeArray)
    oedgesPolyData.GetCellData().SetScalars(oedgeColors)

    tubes = vtk.vtkTubeFilter()
    tubes.SetInput(cutedgesPolyData)
    tubes.SetRadius(5*_rtube)
    tubes.CappingOff()
    tubes.SetNumberOfSides(12)
    tubes.SetGenerateTCoordsToOff()

    otubes = vtk.vtkTubeFilter()
    otubes.SetInput(oedgesPolyData)
    otubes.SetRadius(_rtube)
    otubes.CappingOff()
    otubes.SetNumberOfSides(12)
    otubes.SetGenerateTCoordsToOff()
        
    spheresMap = vtk.vtkPolyDataMapper()
    spheresMap.SetInputConnection(spheres.GetOutputPort())
    spheresMap.SetLookupTable(lut)
    spheresMap.SetScalarRange(1, len(_domains))

    cutedgesMap = vtk.vtkPolyDataMapper()
    cutedgesMap.SetInputConnection(tubes.GetOutputPort())
    cutedgesMap.SetLookupTable(lutEdges)
    cutedgesMap.SetScalarRange(1,1)     # All scalar for cutedges are > to 1, so are going to take the max color of lutEdges (red)
    cutedgesMap.ScalarVisibilityOn()
    oedgesMap = vtk.vtkPolyDataMapper()
    oedgesMap.SetInputConnection(otubes.GetOutputPort())
    oedgesMap.SetLookupTable(lutEdges)
    oedgesMap.SetScalarRange(0, 0)
    oedgesMap.ScalarVisibilityOn()
    
    spheresActor = vtk.vtkActor()
    spheresActor.SetMapper(spheresMap)
    spheresActor.GetProperty().SetSpecularColor(0,1,0)

    cutedgesActor = vtk.vtkActor()
    cutedgesActor.SetMapper(cutedgesMap)
    cutedgesActor.GetProperty().SetOpacity(1.0)
    cutedgesActor.GetProperty().SetInterpolationToGouraud()
    oedgesActor = vtk.vtkActor()
    oedgesActor.SetMapper(oedgesMap)
    oedgesActor.GetProperty().SetOpacity(1.0)
    oedgesActor.GetProperty().SetInterpolationToGouraud()

    ren = vtk.vtkRenderer()
    ren.AddActor(spheresActor)
    if len(_domains)>1:
        ren.AddActor(cutedgesActor)
    ren.AddActor(oedgesActor)

    if fl_legend != 0:
        # Add legend for the domain colors
        for part in _domains:
            text = vtk.vtkVectorText()
            line = "proc "+str(part-1)
            text.SetText(line)
            mapText=vtk.vtkPolyDataMapper()
            mapText.SetInputConnection(text.GetOutputPort())
            textActor=vtk.vtkFollower()
            textActor.SetMapper(mapText)
            textActor.SetScale(scale)
            textActor.AddPosition(xmax+(xmax-xmin)/10+(xmax-xmin)/20,  ymin+(part-1+2)*(ymax-ymin)/10, 0.0)
            if (fl_file == 3):
                textActor.GetProperty().SetColor(0.0, 0.0, 0.0)
            ren.AddActor2D(textActor)
            ren.ResetCameraClippingRange()
            textActor.SetCamera(ren.GetActiveCamera())

        if len(_domains)>1:
            # Legend for the cutting edges 
            text = vtk.vtkVectorText()
            line = "cut e. "#(" + str(len(cutedges)) +")"
            text.SetText(line)
            mapText=vtk.vtkPolyDataMapper()
            mapText.SetInputConnection(text.GetOutputPort())
            textActor=vtk.vtkFollower()
            textActor.SetMapper(mapText)
            textActor.SetScale(scale)
            textActor.AddPosition(xmax+(xmax-xmin)/10+(xmax-xmin)/20, ymin, 0.0)
            if (fl_file == 3):
                textActor.GetProperty().SetColor(0.0, 0.0, 0.0)
            ren.AddActor2D(textActor)
            ren.ResetCameraClippingRange()
            textActor.SetCamera(ren.GetActiveCamera())

        # Legend for the other edges
        text = vtk.vtkVectorText()
        line = "other e. "#(" + str(len(oedges)) + ")"
        text.SetText(line)
        mapText=vtk.vtkPolyDataMapper()
        mapText.SetInputConnection(text.GetOutputPort())
        textActor=vtk.vtkFollower()
        textActor.SetMapper(mapText)
        textActor.SetScale(scale)
        textActor.AddPosition(xmax+(xmax-xmin)/10+(xmax-xmin)/20, ymin+(ymax-ymin)/10, 0.0)
        if (fl_file == 3):
            textActor.GetProperty().SetColor(0.0, 0.0, 0.0)
        ren.AddActor2D(textActor)
        ren.ResetCameraClippingRange()
        textActor.SetCamera(ren.GetActiveCamera())


    if fl_file != 3:
	ren.SetBackground(0.0, 0.0, 0.0)
    else:
	ren.SetBackground(1.0, 1.0, 1.0)

    for label in labels:
        ren.AddActor(label)

    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(ren)
    renWin.SetSize(wx, wy)

    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    style = vtk.vtkInteractorStyleTrackballCamera()
    iren.SetInteractorStyle(style)
    iren.Initialize()


    if (fl_file == 0):                  # Interactive 3D vizualization
        iren.Start()
    elif (fl_file == 1):                # pdf file
        print "Starting normal pdf file generation ..."
        exp = vtk.vtkGL2PSExporter()
        exp.SetRenderWindow(renWin)
        exp.SetFileFormatToPDF()
        #exp.CompressOff()
        #exp.PS3ShadingOff()
        exp.SetFilePrefix(filename)
        exp.Write()
        print filename + ".pdf written!"

    elif (fl_file==2):                  # Raster pdf
        print "Starting raster pdf file generation ..."
        #renWin.Render()
        exp = vtk.vtkGL2PSExporter()
        exp.SetRenderWindow(renWin)
        exp.SetFileFormatToPDF()
        exp.SetFilePrefix(filename)
        exp.Write3DPropsAsRasterImageOn()
        exp.Write()
        print filename + ".pdf written!"
                          
    elif (fl_file == 3):                # png file
        print "Starting png file generation ..."
        renWin.Render()
        w2i = vtk.vtkWindowToImageFilter()
        w2i.SetInput(renWin)
        
        writer = vtk.vtkPNGWriter()
        writer.SetInput(w2i.GetOutput())
        writer.SetFileName(filename + ".png")
        writer.Write()
        print filename + ".png written!"

    elif (fl_file == 4):                # vrml file
        print "Starting VRML file generation ..."
        #renWin.Render()
        writer = vtk.vtkVRMLExporter()
        writer.SetRenderWindow(renWin)
        writer.SetFileName(filename + ".wrl")
        writer.Write()
        print filename + ".wrl written!"


    elif (fl_file == 5):                # obj (Wavefront) file
        print "Starting OBJ file generation ..."
        #renWin.Render()
        writer = vtk.vtkOBJExporter()
        writer.SetRenderWindow(renWin)
        writer.SetFilePrefix(filename)
        writer.Write()
        print filename + ".obj written!"


    elif fl_file == 7:                  # VTK (PolyData) file
        print "Starting VTK file generation"
        writer = vtk.vtkPolyDataWriter()
        writer.SetInput(nodesPoly) 
        writer.SetFileName(filename+".vtp")
        writer.Write()
        print filename + ".vtp written!"
        
    else:
        print "ERROR: Bad file flag. Exit!"
        sys.exit(1)
        

def window(nodes):
    """ Return the window coordinates:
    in order, xmin, xmax, ymin, ymax
    """

    xmin = 100000000
    xmax = -100000000
    ymin = 100000000
    ymax = -100000000
    for node in nodes:
        if node[0]<xmin:
            xmin = node[0]
        elif node[0]>xmax:
            xmax = node[0]
        if node[1]<ymin:
            ymin = node[1]
        elif node[1]>ymax:
            ymax = node[1]

    return xmin-1, xmax+1, ymin-1, ymax+1
        
        

def main():
    """ Main function
    """
    
    gnuplot = "/opt/local/bin/gnuplot"
    temp = 1                            # Delete gnuplot data and script files?
    debug = 0                           # Debug flag
    fl_part_format = 1  # Format of the part file: 0=triangle, 1=metis
    fl_node_format = 0                  # Format of the node and edge file: 0=triangle, 1=gts, 2=gmsh, 3=graph(connectivity graph, nodal or dual)
    fl_render = 0                       # Render option: 0=gnuplot, 1=vtk
    fl_3D = 0                           # Explicitely ask for 3D handling: 1=3D, other values: let the soft decide
    fl_part = 1                         # Draw the partition? 0=no, 1=yes

    fl_paper = 0                        # Draw output adapted for paper: 0=no, 1=yes
    rsphere = 100                       # Size of a sphere representing a node, ptsize for gnuplot
    rtube = 100
    fl_legend = 1                       # Add legend
    label = 0                           # 0=dont print labels, 1=print them

    fl_file = 100                       # Output file type
                                        # VTK type of file:
                                        # 0=No file (3D interactive vizualization) (DEFAULT),
                                        # 1=normal 2D pdf (view with adobe reader),
                                        # 2=raster pdf,
                                        # 3=png
                                        # 4=VRML
                                        # 5=OBJ
                                        # 7=VTK polydata
                                        # Gnuplot type of file
                                        # 1=pdf (DEFAULT)
                                        # 3=png
                                        # 6=svg
    

    in_node = ""
    in_edge = ""
    in_part = ""

    out_gpd = "test.gpd"
    out_gps = "test.gps"
    filename = "test"                   # Filename of the resulting mesh partitioning
    
    try:
        opts, arg = getopt.getopt(sys.argv[1:], "plfo:dhi:n:r:", ["label", "help","debug", "in_part_format=", "in_node_format=", "render", "out_data=", "out_script=", "output=", "paper", "rsphere=", "rtube=", "interactive", "png", "svg", "pdf", "vrml", "obj", "vtk", "rasterpdf", "nolegend", "3D", "nopart"])
    except getopt.error, msg:
        print msg
        print "For help, use --help"
        sys.exit(2)

    # Process options
    for o, a in opts:

        if o in ("-h", "--help"):
            print __doc__
            sys.exit(0)
        elif o == "-t":
            temp = 1
        elif o=="-p":
            temp = 0
        elif o=="--paper":
            fl_paper = 1
        elif o == "--out_data":
            out_gpd = a
            fl_render = 0
        elif o == "--out_script":
            out_gps = a
            fl_render = 0
        elif o in ("-o", "--output"):
            filename = a            
        elif o in ("-l", "--label"):
            label = 1            
        elif o in ("-d", "--debug"):
            debug = 1
            label = 1
        elif o=="--nopart":
            fl_part = 0
        elif o == "--rsphere":
            rsphere = float(a)
        elif o == "--rtube":
            rtube = float(a)
            fl_rtube = 1
        elif o == "--interactive":
            fl_file = 0
            if fl_render == 0:
                print "No interactive visualization option for gnuplot"
                sys.exit(1)
        elif o == "--png":
            fl_file = 3
        elif o == "--pdf":
            fl_file = 1
        elif o == "--rasterpdf":
            fl_file = 2
            if fl_render == 0:
                print "No raster pdf option for render gnuplot"
                sys.exit(1)
        elif o == "--vrml":
            fl_file=4
            if fl_render == 0:
                print "No VRML option for render gnuplot"
                sys.exit(1)
        elif o == "--svg":
            fl_file=6
            if fl_render == 1:
                print "No SVG option for render VTK"
                sys.exit(1)
        elif o == "--obj":
            fl_file=5
            if fl_render == 0:
                print "No OBJ option for render gnuplot"
                sys.exit(1)
        elif o == "--vtk":
            fl_file=7
            if fl_render == 0:
                print "No vtk option for render gnuplot"
                sys.exit(1)
        elif o == "--nolegend":
            fl_legend = 0
        elif o in ("-i", "--in_part_format"):
            if a=="triangle":
                fl_part_format = 0
            elif a=="metis":
                fl_part_format = 1
            else:
                print "Input part format \""+a+"\" unknown!"
                sys.exit(1)
        elif o in ("-n", "--in_node_format"):
            if a=="triangle":
                fl_node_format=0
            elif a=="gts":
                fl_node_format=1
            elif a=="gmsh":
                fl_node_format=2
            elif a=="graph":
                fl_node_format=3
            else:
                print "Input node format\""+a+"\" unknown!"
                sys.exit(1);
        elif o in ("-r","--render"):
            if a=="gnuplot":
                fl_render = 0
            elif a=="vtk":
                fl_render = 1
            else:
                print "Ouput renderer \""+a+"\" unknown!"
                sys.exit(1);
        elif o=="--3D":
            fl_3D = 1
        else:
            print "Option " + o + " unknown. bye ;)"
            sys.exit(2)


    # Default values
    if fl_file == 100:
        if fl_render == 0:
            fl_file = 1                 # Default gnuplot = pdf
        elif fl_render == 1:
            fl_file = 0                 # Default VTK = interactive

    if rsphere == 100:
        if fl_render == 0 and fl_file in (1,3): # pdf or png 
            rsphere = 1
        if fl_render == 0 and fl_file == 6: # svg
            rsphere = 0.6
        if fl_render == 1:
            rsphere = 0.6

    if rtube == 100:
        if fl_render == 1:
            rtube = rsphere/10
        if fl_render == 0:
            if fl_file in (1,3):
                rtube = 1
            elif fl_file == 6:
                rtube = 0.6

        
    # Check input files configuration
    if fl_node_format in (0, 3): # Triangle, nodal or dual format, expect 2 (.node & .edge) or 3 files (.part)
        if len(arg)==2:       # 2 files: must be .node and .edge files
            in_node = arg[0]
            in_edge = arg[1]
            no_draw_part = 1
        elif len(arg)==3:
            in_node = arg[0]
            in_edge = arg[1]
            in_part = arg[2]
        else:
            print "Please give 2 (triangle .node and .ele) or 3 files (including .part file)"
            sys.exit(2);
            
    elif fl_node_format in (1, 2):         # gts format
        if len(arg)==1:
            in_node = arg[0]
            no_draw_part = 1
        elif len(arg)==2:
            in_node = arg[0]
            in_part = arg[1]
        else:
            print "Please give 1 (gts .node pr gmsh .msh) or 2 files (including .part file)"
            sys.exit(2);

    if in_node != "":
        try:
            fnode = file(in_node, "r")
            print "Use "+in_node+" node file"
        except IOError:
            print "Impossible to open the .node file ".in_node
            sys.exit(2)

    if in_edge != "":
        try:
            fedge = file(in_edge, "r")
            print "Use "+in_edge+" edge file"
        except IOError:
            print "Impossible to open the .edge file ".in_edge
            sys.exit(2)

    if in_part != "":
        try:
            fpart = file(in_part, "r")
            print "Use "+in_part+" node file"
        except IOError:
            print "Impossible to open the .part file " + in_part
            sys.exit(2)
    else:
        fl_part = 0
        

    names = [in_node, in_part]

    # Build standard objects
    
    if fl_node_format == 0:
        nodes, dim = pnode(fnode)
        edges = pedges(fedge)
        
    elif fl_node_format == 3:
        nodes, dim = pnode(fnode)
        edges = gedges(fedge)
        
    elif fl_node_format == 1:
        nodes, edges = gts_node(fnode)
        dim = 3
        
    elif fl_node_format == 2:
        nodes, edges = gmsh_node(fnode)
        dim = 3

    else:
        print "Bad node format (must be 0 or 1, is "+fl_node_format+")"
        sys.exit(2)
        
    if fl_part == 1:
        mappings, domains = pparts(fpart, fl_part_format)
    else:
        domains=[1]
        mappings=[(x+1,1) for x in range(len(nodes))]

    if fl_part == 1:
        fpart.close()
    if fl_node_format in (0,3):
        fedge.close()
    fnode.close()

    if debug == 1:
        print "nodes:"
        pprint.pprint(nodes)
        print "edges"
        pprint.pprint(edges)
        print "mappings"
        pprint.pprint(mappings)
        print "domains"
        pprint.pprint(domains)
        
    # Gnuplot rendering
    if fl_render == 0:
        ogpd = file(out_gpd, "w")
        ogps = file(out_gps, "w")
        cutedges = build_gpd(nodes, dim, edges, domains, mappings, ogpd)

        if fl_file in (1,3):                 # pdf output (Normal, paper, ...)
            filename_ps = filename + ".ps"
            build_gps(nodes, dim, edges, domains, mappings, filename, out_gpd, ogps, fl_paper, rsphere, len(cutedges), rtube, fl_file, fl_legend, names)
            ogpd.close(), ogps.close()
            order = "gnuplot %s;\n" % out_gps
            if fl_file == 1:
                filename_pdf = filename + ".pdf"
                order += "ps2pdf %s; rm -f %s" % (filename_ps, filename_ps)
                os.system(order)
                print "%s successfully created!" % filename_pdf
            else:
                filename_png = filename + ".png"
                order += "convert %s %s; rm -f %s" % (filename_ps, filename_png, filename_ps)
                os.system(order)
                print "%s successfully created!" % filename_png
                
            
        elif fl_file == 6:
            filename_svg = filename + ".svg"
            build_gps(nodes, dim, edges, domains, mappings, filename, out_gpd, ogps, fl_paper, rsphere, len(cutedges), rtube, fl_file, fl_legend, names)
            ogpd.close(), ogps.close()
            order = "gnuplot %s;\n" % out_gps
            os.system(order)
            print "%s successfully created!" % filename_svg
            
        else:
            print "Unknown format fl_type=%d" % fl_file
            sys.exit(1)

        if temp == 1:
            os.remove(out_gpd)
            os.remove(out_gps)

    # VTK rendering
    elif fl_render == 1:
        build_vtk(nodes, dim, edges, domains, mappings, rsphere, rtube, filename, fl_file, fl_legend, label)

if __name__ == '__main__':
    main()
    

                                 
    

